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Abstract 

A robust and general solver for the radial Schrodinger, Dirac, and Kohn-Sham equa- 
tions is presented. The formulation admits general potentials and meshes: uniform, 
exponential, or other defined by nodal distribution and derivative functions. For 
a given mesh type, convergence can be controlled systematically by increasing the 
number of grid points. Radial integrations are carried out using a combination of 
asymptotic forms, Runge-Kutta, and implicit Adams methods. Eigenfunctions are 
determined by a combination of bisection and perturbation methods for robustness 
and speed. A novel outward Poisson integration is employed to increase accuracy 
in the core region, allowing absolute accuracies of 10 -8 Hartree to be attained 
for total energies of heavy atoms such as uranium. Detailed convergence studies 
are presented and computational parameters are provided to achieve accuracies 
commonly required in practice. Comparisons to analytic and current-benchmark 
density-functional results for atomic number Z = 1-92 are presented, verifying and 
providing a refinement to current benchmarks. An efficient, modular Fortran 95 im- 
plementation, dftatom, is provided as open source, including examples, tests, and 
wrappers for interface to other languages; wherein particular emphasis is placed on 
the independence (no global variables), reusability, and generality of the individual 
routines. 
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1 Introduction 



Atomic structure calculations are a lynchpin of modern materials theory. They 
provide the basis for understanding the properties of individual atoms and 
play a central role in electronic structure calculations of larger, multi-atom 
systems such as molecules, nanostructures, and solids [1]. In the latter con- 
text, being a core computational kernel for larger-scale calculations, both 
accuracy and efficiency are of prime importance. Such atomic calculations, 
solving radial Schrodinger or Dirac equations, can be found at the heart of 
plasma physics calculations [2] , all-electron electronic structure methods such 
as Korringa-Kohn-Rostoker (KKR) [3], linearized muffin tin orbital (LMTO) 
[4], and linearized augmented planewave (LAPW) [5], pseudo-atomic-orbital 
based methods [6], ab initio pseudopotential construction [1], and projector 
augmented wave (PAW) [7] , relaxed-core PAW [8] , and all-electron pseudopo- 
tential (AEPP) [9] methods, among others. 

Due to their central importance in the full range of electronic structure cal- 
culations, from isolated atoms to solids, a number of atomic structure codes 
have been developed over the decades, see, e.g., [10-15]. However, most have 
been developed within, and remain closely tied to, the specific larger-scale 
method of which they are a part, and hence are difficult to separate and use for 
other purposes. Indeed, it was precisely our own need for a robust and general 
solver for use in a finite-element [16] based AEPP [9] method that motivated 
the present work. Beyond difficulty to extract, however, existing codes have 
widely differing capabilities, efficiencies, robustness, and availability. Table 1 
lists features of several established codes. Some solve the equations of den- 
sity functional theory (DFT) while others solve Hartree-Fock (HF). Some are 
relativistic, others not. Some can be converged to high precision while others 
can be difficult to converge beyond the accuracies required by the larger codes 
of which they are a part. Some allow any mesh while others implement only 
exponential meshes. Some allow any potential while others implement only 
singular all-electron potentials. Some implement perturbation corrections for 
speed, others only bisection. Some are open source, others not. Note that the 
table shows only shooting-type solvers [17], as we discuss in the present work. 
Other non-shooting type approaches (see, e.g., [18,19] and references therein) 
have been developed and implemented as well but have not yet found wide 
adoption in larger-scale electronic structure calculations due in part to robust- 
ness issues (e.g., spurious states) which can arise [19]. Note also that the table 
is intended only as a general guide: in the context of any given feature, "No" 
should be understood to mean only that we did not find it straightforward to 
implement in the referenced code. 

In the present work, we bring together ideas from a host of atomic struc- 
ture codes developed over many decades, and add new ideas to increase ac- 
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Code 


dftatom 


Desclaux 


atompp 


grasp2k 


atompaw 


PEtot 


Elk 


Reference 


This work 


[10] 


[11] 


[12] 


[13] 


[14] 


[15] 


Open Source 


MIT 


No 


No 


No 


GPL 


BSD 


GPL 


Schrodinger 


Yes 


No 


Yes 


No 


Yes 


Yes 


No 


Dirac 


Yes 


Yes 


No 


Yes 


No 


Yes 


Yes 


DFT 


Yes 


No 


Yes 


No 


Yes 


Yes 


Yes 


Any potential 


Yes 


No 


Yes 


No 


No 


No 


Yes 


Any mesh 


Yes 


No 


No 


No 


No 


No 


Yes 


Outward Poisson 


Yes 


No 


No 


No 


No 


No 


No 


Perturb, correct. 


Yes 


Yes 


Yes 


No 


Yes 


Yes 


No 


Bisection fallback 


Yes 


No 


No 


No 


No 


No 


No 



Table 1 

Features of various atomic structure solvers. 

curacy and robustness. We try to do this in such a way as to provide as 
general, flexible, and efficient a tool as possible, facilitating subsequent use 
in a range of applications areas. We develop and implement solvers for the 
radial Schrodinger, Dirac, and Kohn-Sham equations. The formulation allows 
general potentials and meshes: singular Coulomb-, finite pseudo-, or other 
potentials; uniform, exponential, or other meshes, as defined by nodal distri- 
butions and derivatives. For a given mesh type, convergence can be controlled 
systematically by increasing the number of grid points. Radial integrations 
are carried out using a combination of asymptotic forms, Runge-Kutta, and 
implicit Adams methods. Eigenfunctions are determined by a combination of 
bisection and perturbation methods for robustness and speed. We employ a 
novel outward Poisson integration to increase accuracy in the core region, al- 
lowing absolute accuracies of 1CF 8 Hartree to be attained for total energies 
of heavy atoms such as uranium. We show detailed convergence studies and 
provide computational parameters to achieve accuracies commonly required in 
practice. We present comparisons to analytic and current-benchmark density- 
functional results [20,21] for atomic number Z = 1-92, verifying and refining 
current benchmarks. 

We provide an efficient, modular Fortran 95 implementation, dftatom, as 
open source, including examples, tests, and wrappers for interface to other 
languages. The code consists of several clearly arranged Fortran 95 modules, 
designed to be easy to apply for other purposes. The integration routines are 
independent of mesh, which can be specified as desired. Routines for common 
meshes are provided. Extensive tests are included to verify the correctness of 
results and serve as examples for a number of cases: Coulombic, self consistent 
DFT, and non-singular potentials; exponential and hyperbolic meshes; atomic 
numbers Z = 1-92; user-specified occupation numbers; and others. The code 
is opensource, available under the terms of the MIT license, from [22]. C and 
Python wrappers are provided so that the code can be easily used from other 
languages as well as interactively from Python. 
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The remainder of the paper is organized as follows. Section 2 describes the 
radial Schrodinger, Dirac, and Poisson equations, associated asymptotics, and 
DFT equations to be solved. Section 3 describes the numerical methods used 
to solve the radial equations and eigenproblem, and meshes employed. Sec- 
tion 4 shows Schrodinger and Dirac results for a Coulombic potential and har- 
monic oscillator potential, relativistic and nonrelativistic DFT calculations of 
uranium, convergence studies determining parameters for specified accuracies, 
and comparison to established benchmarks [20]. We conclude in Section 5. 



2 Equations 



2.1 Radial Schrodinger equation 



The 3D one-electron Schrodinger equation is given by 

(-|V 2 + V(x))^(x) = ^(x), (1) 

in Hartree atomic units, as we use throughout. For a spherically symmetric 
potential 

V(x) = V(r), (2) 
eigenstates of energy and angular momentum can be written in the form 

^nim(x) = R„l(r) Yl m (^j , (3) 

where n is the principle quantum number, I is the orbital angular momentum 
quantum number, and m is the magnetic quantum number; whereupon it 
follows that R n i(r) satisfies the radial Schrodinger equation 



\r 2 Ki{r)) + (r 2 V + + 1)) R nl (r) = Er 2 R nl (r). (4) 
The functions ipnimi*) and R n i(r) are normalized as 

J |Vw(x)| 2 d 3 z = l, (5) 

'DC 

R 2 nl ( r ydr = l. (6) 



The substitution P n i(r) = rR nt (r) and Q n i(r) = P' n i( r ) = Rniij) + r R' nl (r) can 
be used to write the second-order radial Schrodinger equation as a coupled set 
of first-order equations: 

Ki(r) = Qm(r), (7) 
Q'mir) = 2 (V(r) -E+ P nl (r). (8) 
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Such a first-order formulation facilitates solving both nonrelativistic Schrodinger 
and relativistic Dirac equations using the same techniques. The normalization 
of P(r) is 

poo 

P»dr = l. (9) 



o 



For a potential V(r) which behaves like —Z/r + Vo near the origin, as in the 
vicinity of the nucleus, the asymptotic behaviors of P n i and Q n i are [23] 

P nl (r)=r l+ \ (10) 
Q nl (r) = (l + l)r l . (11) 

For large r, assuming V(r) — > as r — > oo, the asymptotic is [23]: 

P n i(r)=e- Xr , (12) 

Q nl (r) = -\P rd (r), (13) 

where 

A = sp2E~i . (14) 

As solutions of a homogeneous system, P ni and Q n i are unique only up to an ar- 
bitrary multiplicative constant. These small- and large-r asymptotics provide 
starting values and derivatives for inward and outward numerical integrations. 



2.2 Radial Dirac equation 



The one-electron radial Dirac equation can be written as 

P'(r) = -*P(r) + (l^XM. + 2c^j Q(r), (15) 

Q'(r) = - ( g ~ c V(r) ) P(r) + ~Q(r), (16) 

where P(r) and Q(r) are related to the usual large g(r) and small f(r) com- 
ponents of the Dirac equation by 

P(r)=rg(r), (17) 
Q(r)=rf(r). (18) 

See for example [24,25] for a discussion and derivation. The energy E does 
not contain the electron rest mass energy, so it can be compared directly 
to the corresponding nonrelativistic energy of the Schrodinger equation. The 
normalization of P(r) and Q(r) is 



L 



oo 

2/^\ i r^2i 



P\r) + Q\r))dr = l. (19) 

o 
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For potentials of the form V(r) = —Z/r + Z\ + 0(r), the asymptotic behavior 
at the origin for Z ^ (Coulombic) is [26,25] 



where 



P(r) = r p , (20) 
QW-r'^, (21) 



l-f. (22) 



c 

For Z = (nonsingular) the asymptotic is [26], for k < 



P(r) = r l+ \ (23) 

w^ ,+2 |rr|' < 24 > 



and for k > 



Q(r) = r l+1 . (26) 

For large r, assuming V(r) — >■ as r — >■ oo, the asymptotic is [23] 

P(r) = e~ Xr , (27) 

Q(r) = -f ^P(r), (28) 

(29) 

where 



A=^-^ = f2 £ -f. (30) 

As in the case of the Schrodinger equation, P nK and are unique only 
up to an arbitrary multiplicative constant, and the above asymptotics provide 
starting values and derivatives for inward and outward numerical integrations. 

The solutions of the Dirac equation are labeled by quantum numbers n and 
k, where k ^ is an integer. Alternatively, they may be labeled by the triplet 
n, I, s, where s — ±1 distinguishes j — I ± \ states, k is then determined by 

n={ fei = i + |(- = +l). (31 ) 

5 I for j = I - \ (s = -1). V 7 
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2.3 Poisson equation 



The 3D Poisson equation for the Hartree potential V H due to electronic density 
n is given by 

VV H (x) = -47m(x). (32) 
For a spherical density n(x) = n(r), this becomes 

^{r 2 V' H )' = V£(r) + JV» = -47m(r), (33) 
where n(r) is the radial particle (number) density, normalized such that 

//"CO 
n(x) d 3 x = y 47rn(r)r 2 dr. (34) 

Substituting (33) into (34) and integrating, we obtain 

lim r 2 VUr) = -Z, (35) 
from which it follows that the asymptotic behavior of V H {r) is 

V n( r ) = r ^ °°- ( 36 ) 

Integrating (36) and requiring — >■ as r — >■ oo then gives the corresponding 
asymptotic behavior for Vff(r): 

V^(r) = -, r -> oo. (37) 
r 

For small r, the asymptotic behavior can be obtained by expanding n(r) about 
r = 0: n(r) = Y^jLo c j r ^ ■ Substituting into Poisson equation (33) gives 

oo 

(^X)' = -4vrE^ +2 - (38) 
Integrating and requiring Vh(0) finite then gives 

oo r j+l 

V^(r) = -4vr ]Tc,— , (39) 

3=0 3+6 

with linear leading term; so that we have 

V' H {r) oc r, r ->■ 0. (40) 

Integrating (39) then gives 

oo r j+2 

V " (r)= - 4 ^(j^m^ +c ' (41) 



with leading constant term C = Vh(0) determined by Coulomb's law: 

POO 

V H (0) = 4vr / rn(r)dr. (42) 
Jo 

Finally, from (40) we have that 

V£(0)=0. (43) 

The above asymptotics provide starting values and derivatives for inward and 
outward numerical integrations. 

2.4 Kohn-Sham equations 

The Kohn-Sham equations consist of the radial Schrodinger or Dirac equations 
above with an effective potential V(r) = V m (r) given by (see, e.g., [1]) 

Vin = V H + V xc + v, (44) 

where V H is the Hartree potential given by the solution of the radial Poisson 
equation (33), V xc is the exchange-correlation potential and v = — — is the 
nuclear potential. 

The total energy is given by 

E [n\ = T s [n] + E H [n\ + E xc [n] + V[n], (45) 

the sum of kinetic energy 

T s [n] = fmeni ~ 4tt J V hl {r)n{r)r 2 dr, (46) 

nl 

Hartree energy 

E H [n] = 2rr J V H {r)n{r)r 2 dr, (47) 
exchange-correlation energy 

E xc [n] = Air J e xc (r; n)n(r)r 2 dr, (48) 

and Coulomb energy 

V[n] = An J v{r)n{r)r 2 dr = — AnZ j n(r)rdr; (49) 

with electronic density in the nonrelativistic case given by 

1 P 2 (r) 

nl 
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where P n i is the radial wavefunction (Eq. (7)) and f n i the associated electronic 
occupation. In the relativistic case, the electronic density is given by 

n(r) = -^E/^ Pl(r)+ 2 gL(r) , (51) 
Air , r A 

nls 

where P n i s and Q n i s are the two components of the Dirac solution (Eqs. (15), (16)) 
and f n i s is the occupation. In both cases above, n(r) is the electronic particle 
density [electrons/volume], everywhere positive, as distinct from the electronic 
charge density p(r) [charge/volume]: p(r) = —n{r) in atomic units. 

These equations are solved self-consistently (see, e.g., [1]): an initial density 
n m and corresponding potential V in are constructed; the Schrodinger or Dirac 
equation is solved to determine wavefunctions R n i or spinor components P 
and Q, respectively; from these a new density n out and potential V^ ut are 
constructed; and from these, a new input density n in and potential V[ n are 
constructed. The process is continued until the difference of V m and V out and/or 
n in and n ont is within a specified tolerance, at which point self- consistency 
is achieved. This fixed point iteration is known as the self- consistent field 
(SCF) iteration. We employ an adaptive linear mixing scheme, with optimized 
weights for each component of the potential to construct new input potentials 
for successive SCF iterations. In order to reduce the number of SCF iterations, 
we use a Thomas-Fermi (TF) approximation [27] for the initial density and 
potential: 

V(r) = -^£2 (52) 
r 

where 

Z cS {r) = Z(l + a^c + /3xe- 7 ^) V 2Q ^, (53) 



x = r(— ) , (54) 

a = 0.7280642371, (55) 
p = -0.5430794693, (56) 
7 = 0.3612163121. (57) 

The corresponding charge density is then 

p(r) = ~(-2V(r))K (58) 

For simplicity, we consider only local-density approximation (LDA) and rela- 
tivistic local-density approximation (RLDA) exchange and correlation poten- 
tials here. More sophisticated functionals can be readily incorporated. Our 
df tatom implementation uses the same parameterization as in NIST bench- 
marks [20]: 

V xc (r;n) = ^-(ne L x ?(n)), (59) 
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where the exchange and correlation energy density e^f can be written as [1] 

^{n) = e L x D {n) + e L c D {n) ) (60) 
with electron gas exchange term [1] 

e^H = -|-(37r 2 n)i (61) 
and Vosko-Wilkes-Nussair (VWN) [28] correlation term 

ld, \ \, ( <r \ , 26 ' ■ ' ^ 



e,(n) « - In 4r + ^ arctan - 
c 1 J 2 \ \Y(y)J Q \2y + b / 



byo 



Y(yo) 



J^ + 2(H 2 ,) ardaD / Q 



, (62) 



in which y = y/r~ S) Y(y) = y 2 + by + c, Q = y/4c - b 2 , y = -0.10498, 
b = 3.72744, c = 12.9352, A = 0.0621814, and 



(63) 



is the Wigner-Seitz radius, which gives the mean distance between electrons. 
In the relativistic (RLDA) case, a correction to the LDA exchange energy 
density and potential is given by MacDonald and Vosko [29]: 

D (n) = e L x D (n)R, (64) 

R=1 -2{ J 2 J ' 

V x RLD (n) = V x LD (n)S, (65) 
31og(/3 + /i) l 

m ~ 2 ' 



where \i = y/1 + (3 2 and (3 



(37r 2 ra)5 _ 4we^ D (n) 
c 3c 



3 Methods of solution 



In this section, we discuss the radial integration methods, meshes, and eigen- 
function isolation methods employed in the Kohn-Sham solution. 
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3. 1 Runge-Kutta and Adams methods 



To allow general, nonuniform meshes, all methods first transform equations 
on a general mesh R(t), l<t<N + l,to equations on a uniform mesh t with 
step size h — 1. If the solution on a general mesh is P(r) and the transformed 
solution on the uniform mesh is u(t) then 

u(t) = P(R(t)) (66) 

(68) 

The methods below require the values u(t) and derivatives on the uniform 
mesh, which we express in terms of the values P(r) and derivatives P'(r) on 
the general mesh, and derivative R'(t) of the function defining the general 
mesh. 

The Runge-Kutta family of methods for the numerical integration of ordi- 
nary differential equations require the values of dependent and independent 
variables in the interior of the element (step). For example, the 4th-order 
Runge-Kutta (RK4) step for an equation of the form y' = f(x,y) is: 

y t+ i = y i + l(k 1 + 2k 2 + 2k 3 + k 4 ), (69) 
6 

ki = f(xi,yi), (70) 

k 2 = f(x i+ i,Vi + \ki), (71) 



2 



h = f(x i+1 _,y t + lk 2 ), (72) 



2 



h = f{xi+i,Vi + fa)- (73) 

Implicit Adams methods, on the other hand, use an extrapolation formula to 
advance the solution to the next grid point and then an interpolation formula 
to correct the value (predictor-corrector) using the grid points only. The 4th- 
order Adams outward extrapolation is given by 

Vi+i =Vi + ^(55*4 - BQyU + 37y' t _ 2 - 9^_ 3 ) (74) 



and interpolation by 



Vi+i = y l + y a {H + i + !9y- - H-i + vU)- (75) 



Since the Schrodinger and Dirac equations are linear (for given effective poten- 
tial in each SCF iteration), the interpolation can be determined analytically, 
thus eliminating predictor-corrector iterations, speeding up the calculation 
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considerably. Both Schrodinger and Dirac equations can be written as 







fp\ 




= c 




K Q'J 







where for the Schrodinger equation, the matrix C is given by 



C{i) 



and for the Dirac equation, 







1 



\J$ + 2(V(i) -E) y 



C(i) 



( * *tm + 2c\ 

R(i) c 



E-V(i) 



m ) 



Then for outward integration, analytic interpolation gives: 



A 
A 

M 
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24' 

1 + \ 2 R' 2 (i + l)detC(i + 1), 
// 



u 2 (i + 1) 



1 
A 

M 



10 + A* (i+ i) r* (<+i) crf+i) 

\o iy v c 2 i(i + i) -Cn(i + i; 



ui(i) + ^(19<(i) - 5«i(i - 1) + ui(i - 2)) 
u 2 (i) + £(19^(i) - 5u' 2 (i - 1) + «' 2 (i - 2)) ; 



(76) 



(77) 



(78) 



(79) 
(80) 

, (81) 
(82) 



See [23] for a detailed discussion. 

For inward integration, the 4th order Adams inward extrapolation is given by 



Vi-i =Vi- ^( 55 ^ - 59 v'i+i + 37 v'i+2 - H+a) 



(83) 



and interpolation by 



Vi-i =Vi- ^(Qy'i-i + I9y- - H+i + j/i+2)> 



(84) 
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in which case, analytic interpolation gives: 



A 
A 

M = 




9 

~24' 

1 + \ 2 R' 2 (i-l)detC(i-l), 




+ \R'(i - I) 



-C 22 (i-1) C 12 (i-1 
\C 21 (i-l) -C u (t-1))) 



M 



ui(i) - ^(19ui(i) - 5«i(i + 1) + ui(i + 2)) 
v ti2(i) - - 5u' 2 (i + 1) + u' 2 (i + 2)) / 



(85) 
(86) 

(87) 



In practice, we use RK4 to compute the first 4 points needed by implicit 
Adams, then switch to the more accurate implicit Adams method to compute 
the remaining points. 



3.2 Mesh 



The methods discussed here can use any desired mesh: uniform, hyperbolic, 
exponential, or other. In practice, however, as we show below, most atomic 
structure codes use some form of exponential (or "logarithmic" ) mesh, as this 
gives an efficient concentration of mesh points in the vicinity of the Coulomb 
singularity, where wavef unctions, densities, and potentials vary most rapidly. 
Given their particular prevalence, we discuss exponential meshes below, and 
provide associated routines in the dftatom distribution. 

Every exponential mesh can be written as a function of exactly four parame- 
ters, r min , r max , a, and N, as follows: 



Ri = a(e^-V -l)+ rmill , (89) 

'"max 

log a 



a = (90 ) 



N- 1' 



(91) 



for % = 1, 2, . . . , N + 1. Note that in the above we have explicitly included 
an arbitrary shift of origin. Without this, three parameters are of course suf- 
ficient to specify an exponential mesh. As we elaborate below, however, the 
flexibility to choose the initial point independent of other mesh parameters 
affords additional efficiency when employing asymptotic expressions for initial 
values. In the above form, all parameters have direct physical meaning: the 
mesh starts at R\ = r min and ends at Rn+i = r max ; as shown below, the pa- 
rameter a is the ratio of rightmost to leftmost element (mesh interval) lengths 
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(determining the mesh gradation), and N is the number of elements in the 
mesh. The ratio q of any two successive elements can be calculated from (89): 



Rn Rfi—1 



e i J = a ^. (92) 



It follows that the mesh gradation a can be expressed using the ratio q as 
a = q N " 1 and thus a is the ratio of rightmost to leftmost element lengths. 

The main advantage of this parametrization is that r min and r max are decoupled 
from the mesh gradation and number of elements and can thus be varied inde- 
pendently. Once r min is determined to provide sufficiently accurate asymptotic 
values, there is still freedom to change the gradation of the mesh by changing 
a. Then, after determining the optimal a, mesh convergence can be achieved 
by simply increasing N. By having the option to optimize a independent of 
both r min and N, one can reduce the number of elements N required for a 
given accuracy, thus speeding up the calculation. 

Exponential meshes in common use are special cases of the above more general 
one. For example the first and second meshes in [20] are given by 

n = r miQ ( J . (93) 



Here, the meaning of i, N, r min , and r max is the same as in (89). This form is a 

JV-l 

special case of (89) with a = (r max /r min ) N . As such, the mesh is determined 
by three parameters r min , r max , and N, and the mesh gradation a depends on 
iV and the fraction r max /r min . Note that this mesh does not allow r min = 
(the mesh gradation would become infinite). For uranium, the first mesh has 
N = 15788, r min = 16Q 1 x92 , and r max = 50; from which it follows that a = 

7.353705 x 10 5 . The second mesh has N = 8000, r min = and r - 800 • 



92 ' ' max v/92 ' 

from which it follows a = 7.651530 x 10 9 . 

The third mesh in [20] is given by 

r t = A(e m - 1). (94) 

The meaning of % is the same as in (89). The mesh is determined by three 
parameters A, B, and N. This is a special case of (89) with r min = A(e B - 1), 
r max = A (e B( - N+1 ^ — 1), and a = e B ( N ~ 1 \ For uranium, the values in [20] are 



A = ^ x 10~ 6 , B = 0.002304, and N = 9019; from which it follows that 
r min = 1.088140 x 10~ 10 , r max = 50.031620, and a = 1.055702 x 10 9 . 

The fourth mesh in [20] is given by: 

Pi = log Ri, (95) 
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^min '"max 



a 



N 



LDA 1.0e-07 50.0 2.7e+06 4866 
RLDA 1.0e-08 50.0 6.2e+07 5268 

Table 2 

Optimal mesh parameters for 10~ 6 a.u. accuracy in total energy of uranium. 

where pi is a uniform mesh with N elements such that R± = r min and Rn+i = 
r max . It follows that pi = ^ log rma * + logr min . Then substituting into (95) 

^miii 

and simplifying, the expression (93) is obtained. As such, the fourth mesh is 
equivalent to the first. For uranium, the values in [20] are N = 2837, r min = 
^2~, and r max = 50; from which it follows that a = 4.564065 x 10 9 . 

To summarize, the parameters of the four meshes for uranium in [20] are then 
as follows. 

I r_min r_max a N 



1 I 6.79347826087e-05 50.0 7.353705e+05 15788 

2 | 1.08695652174e-08 83.4057656228 7.651530e+09 8000 

3 I 1.08814001246e-10 50.0316203306 1.055702e+09 9019 

4 | 1.08695652174e-08 50.0 4.564065e+09 2837 

In Section 4.3, we determine optimal parameters for LDA and RLDA calcu- 
lations, as shown in Table 2. Here, the given r min , r max , and a can be used 
for all atoms, while the given N is sufficient for 10~ 6 a.u. accuracy in total 
energy for uranium; it can be increased for higher accuracy and decreased for 
lower atomic numbers Z to improve speed while retaining 10~ 6 accuracy, as 
we show in Section 4.3. 

3.3 Shooting method for Schrodinger & Dirac eigenproblems 

The required eigenfunctions are determined by the shooting method [17]: i.e., 
eigenvalues are guessed, radial integrations are performed, guesses are up- 
dated, and the process is repeated until convergence is achieved to the desired 
tolerance. For efficiency, we employ a combination of inward and outward inte- 
grations, bisection for robustness, and perturbation theory for efficiency once 
the solution is sufficiently close to convergence. 

The outward integration starts at r min using the asymptotic from Section 2 
with r = r mm , then RK4 is used for the first four steps, then Adams for the 
rest. For the inward integration, a starting point r m is determined such that 
e -Ar m j g on ^ ne orc [ er Q f machine epsilon (e.g., ~ 10~ 16 in double precision), 
with A determined from the large-r asymptotic for the Schrodinger or Dirac 
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equation, as appropriate (Section 2). Then e Arm is used for the last four points 
up to r m and inward Adams is used for the rest. 

The outward integration is performed up to the classical turning point r ctp , 
defined by V(r ctp ) = E. Bisection is used to converge the energy to the point 
that the associated radial solution has the correct number of nodes (n — l — 1). 
After the correct number nodes is obtained, two approaches are considered to 
complete the eigenfunction computation. 

The first approach is bisection: the outward integration is continued from r ctp 
to the rest of the domain and the total number of nodes is counted. If the 
number is strictly greater than n — l — 1, then the energy E is above the 
eigenvalue, otherwise it is below the eigenvalue. Bisection is used to refine the 
eigenvalue to specified accuracy. 

The second approach uses a perturbation correction for the energy, using both 
inward and outward integrations. It typically converges an order of magni- 
tude faster than outward-only bisection once the correct number of nodes is 
obtained. For the Schrodinger equation, the following correction is used [23]: 

E *~ E * + 2/~i*(r) dr ' (%) 

where E 1 is the energy used to calculate P(r) and Q(r). The inward integration 
is multiplied by a constant so that it matches (in value) the outward integra- 
tion at the classical turning point r ctp . This leaves a jump in the derivative of 
the outward integration Q(r~ tp ) and inward integration Q(r^ tp ). E2 is the new 
energy for the next iteration. Usually around five iterations are sufficient to 
achieve 10~ 13 accuracy in energy. 

For Dirac equation, the following correction is used [23]: 
F ~ F 1 /(^ P )(Q(rctp)-g(r c + tp )) 



After the correct number of nodes is determined, the perturbation correction is 
tried first. If it fails to converge (as can occur, e.g., for a too-confined domain), 
we fall back to outward-only bisection for robustness. For positive energies, 
the perturbation correction is inapplicable and outward-only bisection is used 
exclusively. Physically, the outward-only bisection method solves the atom in 
an infinite potential well of the size of the domain while the inward-outward 
perturbation method solves an atom in an infinite domain (by virtue of the 
asymptotic inward integration). As such, for large r max , the two methods con- 
verge to the same value, but for small r max the eigenvalues can differ between 
methods. (An r max study is needed to determine sufficient r max , as we show in 
Section 4). 
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All normalization and other integrals are calculated on [r min ,r max ], avoiding 
the Coulomb singularity at r = 0. Hence, convergence with respect to r min 
(typically ~ lCT 8 a.u.) must be checked, as we show in Section 4. 

3.4 Outward Poisson integration 

Since asymptotics are known at both ends of the domain (Section 2.3), it 
is possible to integrate the Poisson equation either inward or outward. We 
employ outward integration here in order to more accurately resolve rapid 
variations in core region. This better resolves the most tightly bound states 
(e.g., Is and 2s), in particular. 

Since charge is omitted around the Coulomb singularity at r = (0 < r < 
?"min), we take as initial conditions, consistent with Eqs. (42) and (43): 



Then RK4 is used for the first four points, after which the Adams method 
is used with predictor-corrector for the rest. Here again, convergence with 
respect to r min must be checked, as we show in Section 4. 

Having a precise Poisson solver is crucial to resolve all Kohn-Sham states 
accurately. With this outward solver, Kohn-Sham eigenvalues are typically 
an order of magnitude more precise than the total energy, with the most 
pronounced effect on the more tightly bound states. This is in contrast to 
other solvers [20] which typically resolve eigenvalues less accurately than total 
energies. 



4 Results 

4-1 Analytic test case: V = —Z/r 

To verify the accuracy of the solver and establish convergence indicators, the 
Schrodinger and Dirac equations are solved for potential V = —Z/r with 
Z = 92, for which analytic results are available. r min , r max , and iV convergence 
studies are run for all eigenvalues with n < 7. In the resulting plots, E — E prcv 
gives the change in energy with change of parameter being studied: e.g., as 
r min is decreased toward convergence, r min = 1CT 4 , 10~ 5 , . . . , 10~ 14 , E — E prev 
at r min = 1CT 7 is the difference of computed energies at r min = 1CT 7 and 




(98) 
(99) 



V^(r min ) = 0. 
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Table 3 

Computed and exact Schrodinger eigenvalues for V = — — . 

r m i n = 10~ 6 . E — E conv is the difference of the computed energy from the fully 
converged value: e.g., E — E conv at r min = 1CT 7 is the difference of energy 
at r min = 1CT 7 from the fully converged value at r min sufficiently small that 
reducing it further does not reduce the error further (due to finite precision). 

For the Schrodinger equation, the r m j n convergence study was run with a = 
10 9 , r max = 50, and iV = 50000. From the resulting plot (Fig. la), we see 
that r min = 10~ 7 is a converged value; decreasing r min further does not reduce 
energy differences further. The r max convergence study was run with a = 
2.7 • 10 6 , r min = 10~ 7 , and N = 50000. From the resulting plot (Fig. lb), 
we see that r max = 3.5 is a converged value. Finally, the N study was run 
with a = 2.7 • 10 6 , r min = 10~ 7 , and r max = 50. From the plot (Fig. lc), it is 
determined that N = 50000 is a converged value. 

From these results, it is seen that the fixed parameters for each convergence 
study were fully converged and the error indicator in each case was below 10~ 10 
a.u. Table 3 compares the computed eigenvalues to exact values (—^2) for this 
analytic test case and it is verified that the absolute error of the computed 
values is below 10~ 10 , as indicated by the convergence studies. 

Of course, smaller iV could be used to achieve the above accuracies using the 
converged r min and r max values indicated by the above studies, then optimizing 
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rmax study for the sum of eigenvalues 













nhstF- 


F. 1 































































































































































































































































■ 1? — 

0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 



max 

(b) r max study 



N study for the sum of eigenvalues 



m 10 a 

10" 9 











-- «>»(E-E W „) 













































































































10000 20000 30000 40000 50000 60000 



N 

(c) N study 

Fig. 1. Convergence studies for the sum of Schrodinger eigenvalues for V = —Z/r. 



the mesh gradation a to minimize the required N. The purpose of the present 
analytic case study, however, is to verify the accuracy of the solver and es- 
tablish the robustness of convergence indicators, E — E prcv and E — E conv . 
Optimal a and iV are discussed in Section 4.3. 

For the Dirac equation, the iV study is first done for r min larger than 10 -11 
and iV = 200000 was obtained as a converged value (a = 10 9 was used so that 
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Table 4 

Computed and exact Dirac eigenvalues for V = — j . 



convergence was obtained for tractable for all r m i n ). The r m i n study was 
then run and determined r min = 10~ 8 as a converged value, giving accuracy 
better than 5 • lCT 10 a.u. r max = 50 was used and it was verified that it is a 
converged value. Table 4 shows a comparison of computed and exact Dirac 
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eigenvalues 



c 2 



E nK = - c\ 100 

ft | (Z/cY 
V («-|«|+/3) 2 

/3 = - (Z/c) 2 , (101) 

for the Coulomb problem, and it is verified that the absolute error is below 
5 • 10~ 10 a.u., as indicated by the convergence studies. Note: all numerical and 
analytic calculations in this paper use the 1986 CODATA [30] recommended 
value c = 137.0359895 a.u., the same value as in NIST benchmarks [20]. 



4-2 Nonsingular test case: V = \uj 2 r 2 



Table 5 shows computed Schrodinger eigenvalues for the harmonic oscillator 

V 



potential V(r) = \oj 2 r 2 with ou — 1 compared to exact values 



E nl =u(2n-l-^. (102) 

Asymptotics (10), (11) are used for outward integration, with perturbation 
correction turned off due to positive energies. It is seen that r m i n = 10~ 7 a.u., 
r max = 10 a.u., a = 10 and N = 5000 are sufficient to obtain eigenvalues 
accurate to < 10~ 8 a.u. 

Table 6 shows Dirac eigenvalues for the harmonic oscillator potential V(r) = 
7jU 2 r 2 with u — 1. Asymptotics (23), (24) and (25), (26) are used for outward 
integration, with perturbation correction again turned off due to positive en- 
ergies. Lacking analytic results for comparison, we compared instead to inde- 
pendent high-order finite-element based calculations [31] which require no as- 
sumptions of asymptotic forms or determination of sufficiently small r min > 0. 
For r m i n = 10~ 8 a.u., r max = 10 a.u., a = 80, and = 5000, we find agreement 
of all eigenvalues to < 10~ 8 a.u., with values very near to the corresponding 
Schrodinger results due to the relative weakness of the potential, in stark con- 
trast to the singular Coulombic case. Note that, unlike the case of a singular 
potential, the number of nodes for k > is only n — I, whereas the number of 
nodes for k < remains n — I — 1, as in the singular case. 



4-3 Uranium 



We now consider full, self-consistent nonrelativistic (LDA) and relativistic 
(RLDA) Kohn-Sham calculations of uranium: a stringent test case, as numer- 
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Table 5 

Computed and exact Schrodinger eigenvalues for V = \r 2 . 
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5.49976206 


7 


1 





12.49847782 


4 


2 


1 


5.49969551 


7 


1 


1 


12.49843790 


4 


3 





4.49989517 


7 


2 





11.49874396 


4 


3 


1 


4.49980199 


7 


2 


1 


11.49867742 


5 








9.49911657 


7 


3 





10.49899680 


5 


1 





8.49931620 


7 


3 


1 


10.49890364 


5 


1 


1 


8.49927627 


7 


4 





9.49923634 


5 


2 





7.49950252 


7 


4 


1 


9.49911657 


5 


2 


1 


7.49943598 


7 


5 





8.49946258 


5 


3 





6.49967554 


7 


5 


1 


8.49931619 


5 


3 


1 


6.49958238 


7 


6 





7.49967553 


5 


4 





5.49983526 


7 


6 


1 


7.49950251 


5 


4 


1 


5.49971547 











Table 6 

Computed Dirac eigenvalues for V = \r 2 . 
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ous eigenstates are required, with tightly bound, highly oscillatory s states, 
spanning energies from ~ —0.1 to ~ —40,000 a.u. 

We first determine sufficiently large r max that the associated error (due to 
confining wavefunctions and potentials) is below the level of other sources in 
the calculation. We start by setting r m j n , r max , a, and N to typical values 
consistent with previous findings [20]. We then increase N to converge to the 
limit of finite precision. Then, increasing r max in LDA (Figs. 2b and 2c) and 
RLDA (Figs. 3b and 3c) calculations (with perturbation correction turned 
off), considering both eigenvalues and total energy, shows that r max = 50 a.u. 
is a well converged value for both LDA and RLDA, consistent with previous 
findings [20]. 

The aim of the next N study is then to find N (and a) such that the total 
energy is converged for broad range of r m i n , so that such N can be used for the 
r m i n study. In theory, by simply increasing N, the total energy must eventually 
converge, but for small a (low concentration of grid points around the origin), 
the required N might be very high. There is some coupling of r min and a due 
to closer approach to the singularity at r = for smaller r min : smaller r min 
requires somewhat larger a for a given N and accuracy. It was found that using 
a = 10 9 , full convergence is achieved for all r min > 10 -14 with N = 50000. 

The initial point r min > determines the accuracy of the asymptotics used to 
start the radial integrations, which are exact only in the limit r — > 0. The goal 
of the r m i n study is to find fully converged r m ; n such that the associated error 
is at the limit of finite precision; so that decreasing further does not improve 
results. To do so, the fully converged a and N are taken from the initial iV 
study, then converged total energies are calculated for all r min from 10~ 14 to 
10" 5 . 

Fig. 2a shows the r min study for the LDA calculation. As can be seen, the 
total energy is converged for r min = 10~ 7 and decreasing r m i n further, the error 
indicators remain below the 5-10~ 9 level. Choosing r min smaller than 10~ 7 does 
not increase accuracy further. In general, choosing r min as large as possible 
avoids the need for excessive grid points to resolve rapidly varying potentials, 
densities, and wavefunctions in the vicinity of the Coulomb singularity at 
r = 0. As such, r min = 10~ 7 is used, as it gives sufficiently accurate asymptotics 
such that the associated errors are reduced to the limits of finite precision 
(which eliminates r min from consideration when choosing other parameters), 
is consistent with values adopted previously [20], and is not so small as to 
require excessive grid points around the origin. 

Fig. 3a shows the r min study for the RLDA calculation. Using the same pro- 
cedure as for the LDA calculation, r min = 10~ 8 was chosen. 

Now the optimal mesh for 10~ 6 a.u. accuracy in total energy is found by 
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rmax study for total energy _ 3 N study for total energy 




N 



(c) r max study total energy (d) N study total energy 

Fig. 2. LDA (Schrodinger) DFT convergence studies for U. 

setting r min and r max to converged values and varying the mesh gradation a to 
minimize the number of elements N required to achieve the specified accuracy. 
Proceeding in this way, we find optimal mesh parameters as in Table 2 for 
1CT 6 a.u. accuracy in the total energy of uranium. 

Having determined converged r min , r max , and optimized a (Table 2), we take 
N to convergence (Fig. 2d) to find converged total energy 

E tot = -25658.41788885 a.u. (103) 

and orbital energies (Table 7) with error < 1CT 8 a.u. for the LDA calculation. 

Proceeding similarly for the RLDA calculation, we take N to convergence 
(Fig. 3d) to find converged total energy 

Etot = -28001. 13232548 a.u. (104) 

and orbital energies (Table 8) with error < 10~ 8 a.u. for the RLDA calculation 
also. 
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state 


eigenvalue [a.u.] 


Is 


-3689 


35513984 


12s 


-639 


77872809 


2p 


-619 


10855018 


3s 


-161 


11807321 


3p 


-150 


97898016 


3d 


-131 


97735828 


4s 


-40 


52808425 


4p 


-35 


85332083 


4d 


-27 


12321230 


4f 


-15 


02746007 


5s 


-8 


82408940 


5p 


-7 


01809220 


5d 


-3 


86617513 


5f 


-0 


36654335 


6s 


-1 


32597632 


6p 


-0 


82253797 


6d 


-0 


14319018 


7s 


-0 


13094786 



Table 7 

Computed LDA eigenvalues for uranium (10 -8 a.u. accurate). 
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state 


eigenvalue [a.u.] 


ls l/2 


-4223.41902045 


2s l/2 


-789.48978233 


2 J>3/2 


-761.37447597 


2 Pl/2 


-622.84809456 


3s l/2 


-199.42980564 


3 J>3/2 


-186.66371312 


3P1/2 


-154.70102667 


3d 5 / 2 


-134.54118029 


3d 3 / 2 


-128.01665738 


4s l/2 


-50.78894806 


4 J>3/2 


-45.03717129 


4 Pl/2 


-36.68861049 


4^5/2 


-27.52930624 


4d 3/2 


-25.98542891 


4/7/2 


-13.88951423 


4/5/2 


-13.48546969 


Bsi/2 


-11.29558710 


5P3/2 


-9.05796425 


5P1/2 


-7.06929564 


5d 5/2 


-3.79741623 


5 d 3/2 


-3.50121718 


5 /7/2 


-0.14678839 


5/5/2 


-0.11604717 




-1.74803995 


6j>3/2 


-1.10111900 


6P1/2 


-0.77578418 


6d 5 / 2 


-0.10304082 


6^3/2 


-0.08480202 


7s l/2 


-0.16094728 



Table 8 

Computed RLDA eigenvalues for uranium (lCT 8 a.u. accurate). 
4.3.1 Optimal N 

In the previous section, for converged r min , r max , and optimized a (Table 2), 
sufficient iV was determined to achieve 10~ 6 accuracy in total energies for 
uranium. In this section, we determine sufficient iV for a series of accuracies, 
from 1CT 3 to 1CT 8 a.u., for all elements Z — 1 - 92. 

Fig. 4a shows the required iV for LDA (Schrodinger) and Fig. 4b for RLDA 
(Dirac) Kohn-Sham equations. 

For simplicity, we use the same converged r min , r max , and optimized a (Table 2) 
for all atoms and just vary for each atom to find a sufficient value for the 
specified accuracy. 

The A^ shown in Figs. 4a and 4b is mainly for use as a starting point. Sufficient 
N will vary somewhat with the choice of exchange-correlation, and very much 
so with the choice of r min , r max , and a. If any of these are changed, as may 
be desired to find more optimal values for different atoms and/or accuracy 
requirements, a new A^ study will be required to determine sufficient N. If, 
however, the converged r min , r max , and a from Table 2 are used, then con- 
vergence can be attained for any atom to the limits of machine precision by 
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(a) Schrodinger (LDA) (b) Dirac (RLDA) 

Fig. 4. Accuracy graph for LDA and RLDA DFT calculations, 
simply increasing N. 

By setting N = 50000 with converged r min , r max , and a from Table 2, fully 
converged total energies and Kohn-Sham eigenvalues are obtained, and used 
to compute errors for Figs. 4a and 4b. 

4-3.2 Comparison to current benchmarks 

Using the meshes from Table 2 with converged N = 50000, the computed 
Kohn-Sham total energies and eigenvalues were compared with current bench- 
mark values [20,21] for all atoms Z = 1 - 92 for LDA (Schrodinger) and RLDA 
(Dirac) Kohn-Sham calculations. It was found that the computed total en- 
ergies were all within 5.29 ■ 10~ 7 of benchmark values for LDA and within 
5.31 ■ 10~ 7 for RLDA; and that the computed eigenvalues were all within 
1.24 x 10~ 6 of benchmark values for LDA and within 1.37 x 10~ 6 for RLDA. 

Reference [20] claims 10~ 6 a.u. absolute accuracy for Kohn-Sham total en- 
ergies and finds maximum deviations among different codes of 2 ■ 10~ 6 for 
associated eigenvalues. Our results, converged to < 10~ 8 a.u., verify that the 
total energies reported in [20] indeed have absolute accuracy 10~ 6 a.u. or bet- 
ter; and furthermore, that the eigenvalues in [20] in fact have absolute accuracy 
2 • 10~ 6 a.u. or better. 



5 Summary and conclusions 

We have presented a robust and general solver for the solution of the radial 
Schrodinger, Dirac, and Kohn-Sham equations of density functional theory; 
and provided a modular, portable, and efficient Fortran 95 implementation, 
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df tatom, along with interfaces to other languages and full suite of examples 
and tests. The solver brings together ideas from many different codes devel- 
oped over the decades, and some new ideas such as outward Poisson integra- 
tion for increased accuracy in the core region and perturbation with fallback 
to bisection for speed and robustness. The solver can accommodate any po- 
tential, whether singular Coulombic or finite, and any mesh, whether linear, 
exponential, or otherwise. We have demonstrated the flexibility and accuracy 
of the associated code with solutions of Schrodinger and Dirac equations for 
Coulombic and harmonic oscillator potentials; and solutions of Kohn-Sham 
and Dirac-Kohn-Sham equations for the challenging case of uranium, obtain- 
ing energies accurate to 1CT 8 a.u., thus verifying and refining by two orders 
of magnitude current benchmarks [20]. We have shown detailed convergence 
studies in each case, providing mesh parameters to facilitate straightforward 
convergence to any desired accuracy by simply increasing the number of mesh 
points. 

At all points in the design of the associated code, we have tried to emphasize 
simplicity and modularity so that the routines provided can be straightfor- 
wardly employed for a range of applications purposes, while retaining high 
efficiency We have made the code available as open source to facilitate dis- 
tribution, modification, and use as needed. We expect the present solvers will 
be of benefit to a range of larger-scale electronic structure methods which rely 
on atomic structure calculations and/or radial integrations more generally as 
key components. 
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